I have made significant changes to BaseBinaryStar::Evolve(). In its current form, even when the option --emit-gravitational-radiation is not selected, COMPAS v03.01.00 creates a different detailed output file file from COMPAS v03.00.05 for the command:

./compas -n 1 --random-seed 0 --logfile-type csv --detailed-output

The value of dt is calculated differently in v03.01.00 even when --emit-gravitational-radiation is not selected. My changes reinstate the code so that v03.01.00 and v03.00.05 calculate dt in the same way (and produce identical detailed output files) when --emit-gravitational-radiation is not selected.

I believe my changes preserve the dt calculation for GW emission, but you should check this. I have checked for the command shown above, and the detailed output files when --emit-gravitational-radiation is selected are identical - EXCEPT that I have changed the code to reduce the duration of the first timestep (see comments in code). If you remove the reduction of the first timestep, the detailed output files are identical my updates vs your original code (if you choose to remove the reduction of the duration of the first timestep, you should only remove it for the GW case - it should be retained for the non-GW case).

As part of those changes I also changed BaseBinaryStar::ChooseTimestep().

I do not like this line in the code:

m_SemiMajorAxis = utils::Compare(aNew, 0.0) > 0 ? aNew : 1E-20;  // if <0, set to arbitrarily small number

Can we not avoid this kludge? Perhaps if we change the definition of HasStarsTouching() from

return (utils::Compare(m_SemiMajorAxis, 0.0) > 0) && (m_SemiMajorAxis <= RSOL_TO_AU * (m_Star1->Radius() + m_Star2->Radius()));

to

return (utils::Compare(m_SemiMajorAxis, 0.0) >= 0) && (m_SemiMajorAxis <= RSOL_TO_AU * (m_Star1->Radius() + m_Star2->Radius())); 

I don't think that will have adverse effects in other places that HasStarsTouching() is used - thoughts?